###############################################################
# Outcomes Mexico 2012 Panel Study design-based approach
###############################################################

rm (list = ls())

library(lmtest)
library(sandwich)
library(stargazer)

# load function that does clustered SEs
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<40){
    warning("Fewer than 40 clusters, variances may be unreliable")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

#########################
# Load and explore data
#########################

# Generate log  
#sink("010_outcomes_mexico_2012_db.txt")

# Read data
load("009_panel_mexico_2012_db_socio_2017jan30.RData")
d_socio = d
d_socio_match = d_match
d = NULL
d_match = NULL

# Read data
load("010_panel_mexico_2012_db_ego_2017jan30.RData")
d_ego = d
d_ego_match = d_match
d = NULL
d_match = NULL

#################################
# Results design based approach
#################################

mod_socio_db = lm(d_socio_match$outcome  ~ d_socio_match$t_ind + 
                                           d_socio_match$vote_incumbent + 
                                           d_socio_match$ideology_imp + 
                                           d_socio_match$vazquez_thermometer_imp +
                                           as.factor(d_socio_match$precint))
mod_db_socio_c <- coeftest(mod_socio_db, vcov = vcovCluster(mod_socio_db, cluster = d_socio_match$precint))

mod_ego_db = lm(d_ego_match$outcome  ~ d_ego_match$t_ind + 
                                       d_ego_match$vote_incumbent + 
                                       d_ego_match$ideology_imp + 
                                       d_ego_match$vazquez_thermometer_imp +
                                       as.factor(d_ego_match$precint))
mod_db_ego_c <- coeftest(mod_ego_db, vcov = vcovCluster(mod_ego_db, cluster = d_ego_match$precint))

##############
# Table LaTeX
##############

# TABLE 5 SECOND PART MANUSCRIPT 
dim(d_socio_match)
dim(d_ego_match)
stargazer(mod_db_socio_c,mod_db_ego_c,
          type = "text",
          title="Design-Based Approach",
          align=TRUE, 
          omit.stat=c("LL","ser","f","rsq","adj.rsq"), 
          no.space=TRUE,  
          multicolumn = TRUE,
          table.placement = "H",
          covariate.labels=c("Negative sociotropic treatment","Negative egotropic treatment"),
          dep.var.caption  = "Voting for the incumbent",
          dep.var.labels.include = FALSE,
          omit = c("vote_incumbent", 
                   "ideology_imp", 
                   "pan_thermometer_imp",
                   "vazquez_thermometer_imp",
                   "precint",
                   "Constant"),
          add.lines = list(c("Controls","Yes","Yes"),
                           c("Neighborhood Fixed effects","Yes","Yes"),
                           c("Observations","60","60")))

#sink()
